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, Abstract: The perfectly matched layers (PMLs), as a boundary termination over an unbounded 

spatial domain, are widely used in numerical simulations of wave propagation problems. Given a 
' set of discretization parameters, a procedure to select the PML profiles based on minimizing the 

. discrete reflectivity is established for frequency domain simulations. We, by extending the function 

(-H I class and adopting a direct search method, improve the former procedure for traveling waves. 
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, 1. Introduction 

o\ . ^ 

\^ , Perfectly matched layers (PMLs) [1] are widely used as boundary terminations in problems 

I to be solved over an unbounded spatial domain. Physically, this is an approach to surround a 

numerical problem domain with a layer of a material which creates as little numerical reflection as 
possible, while also attenuating waves that enter from the problem interior. Mathematically, when 
I a PML is used to truncate the x axis, x is actually replaced by x = /'^ (1 + ia{T))dT, where o" is a 

' real function satisfying certain conditions. 

Nevertheless, in actual numerical simulations, due to the flnite thickness of the PML, reflections 
^ , occur when plane waves incidence upon the PML. Note that PML is generally problem-dependent, 

^ ' for example, the reflection is dependent on the discretization scheme. For frequency domain sim- 

ulations, based on minimizing the average discrete reflectivity, Ya Yan Lu [2] gave a practical 
procedure for selecting the optimal PML proflle a where the function class is restricted to simple 
powers. In this paper, to further reduce the average reflectivity for traveling waves, we consider 
a rational function class for a; to avoid time-consuming computations, we use the Nelder-Mead 
simplex method to determine the coefficients of a. Besides, we also give a simpler a which is easy 
to optimize. Numerical simulations demonstrate that our improved procedure is much better and 
practical. 



2. Motivations to Improve the Former Procedure 

For traveling waves, we consider a two-dimensional waveguiding structure in [2]. The y— component 
of the electric field of a transverse wave satisfies a Helmholtz equation. Consider that the structure 
is unbounded in the negative x direction and the medium is homogeneous (refractive index n = uq) 



for X < G, we then truncate the negative x axis by a PML. For D < H < G, we define (j(x) such 
that cr{x) = for x > H and it(x) > for x < H while cr{H) = and cr'{H) = 0. Replacing x by 
X = (1 + ia{t))dt gives rise to 

s-'^d:,{s-'^d^u) +d'^^u + klnlu = (1) 

where /cq is the free space wavenumber, and the time dependence is e"*"^*. At x = D, we use a 
simple zero boundary condition: u = 0. The actual PML is the layer D < x < H . For H < x < G, 
(1) has a plane wave solution 

u = e^(-"^+'3^) + i?e^("^+^^) (2) 

where the second term is the reflected wave due to the incidence of plane waves upon the PML 
with reflection coefficient R. When x is discretized, R depends on o" [3]. By a second-order finite 
difference approximation in the transverse direction x [2], we can easily find R exactly by solving 
a linear equation system. 

Thus, to select the optimal PML profile a is to find a a such that the following \R\ is minimized: 



— 2 r/"^ 

\R\ = - \R{e)\de. (3) 



For convenience, let 



vr Jo 



= (4) 



In Lu's procedure [2], since a is limited to be a simple power, i.e. 

a = StP, (5) 

we only need to determine p and the dimensionless scaling parameter S such that is minimized. 
The optimal values of S were computed for p = 2, ...,5 respectively, among which the one giving 
the least was chosen to be the overall optimal PML profile. 

In fact, the numerical result can be more satisfying if we give up the restriction of a to be 
simple powers. In practice, people usually use a = for the PML profile [3]. This gives us a 
start to find a better profile. 

After numerical computation and a little adjustment, a = turns out to be better in this 
situation. So in this paper, we first investigate into this rational function class: 



a = —,ap>^,p>2 (6) 



where its coefficients are to be optimized to give a minimal value of the average discrete reflectivity 
Due to the condition a > 0, we let Op > for simplicity while losing certain generality. 
However, if the PML profile is defined by such rational function, the work to determine the 
optimal values of its coefficients becomes much more time-consuming. To save time, we adopt the 
Nelder-Mead (NM) simplex method due to the fact that the problem is nonlinear and the derivative 
information of \R\ is unavailable. 



It is natural to ask why we choose this method. We offer two answers. First, in the NM method, 
the objective function is evaluated at the vertices of a simplex, and movement is away from the 
poorest value. This method tends to work so well in practice by producing a rapid initial decrease 
in function values. Second, when we consider a simpler a with two parameters, the NM algorithm 
gives answer in a much shorter time. Thus, due to its powerful local descent property, we decide 
to adopt this method though it may not give globally optimized solution. ^ 

Note that the NM method is for unconstrained problems, but > in (6) is actually a 
constrainment. To overcome this, we define a as below 

|a2|r^ + laalr^ + ... + |ap|rP 
a = — ,ap >0,p>2 

1 + T 

When similar situation occurs, we will tackle it in this way again without further remarks. 

Furthermore, after optimization for (6), we try to define o" by a simpler function class which 
largely conserves the good properties of the former one but is much easier to optimize and more 
useful in practice. 



(7) 



3. The Improved Procedure with Numerical Results 

A. General Results 

In numerical simulations, we adopt the conditions in [2] so that we could compare the average 
discrete reflectivity between the two procedures. Explicitly we have the wavelength Xq = 1 fi 
m, ko = 27r/Ao, no = 1. 

Consider an example in [2]: A PML with thickness of five grids (m = 5), where the grid size 
h = ^Ao = 0.05 H m. Under the restriction on a in (5): a = St^,p < 5, the author obtained an 

optimal profile: p = 3 and S = 100.4, which gave rise to = 0.013. 

Through our improved procedure, we obtain a simple and better result: a = (23. 6t^ + 
35.9r^)/(l — r) which gives rise to \R\ = 0.0047, only 36% of the former one. Moreover, if we 
give up the simplicity and allow a higher order of the rational function, we could further get an 
result of 0.0031, only 24% of the former one. 

B. A Rational Function Class for the PML Profile 

First of all, we investigate into (6). Take note that the conditions cj(0) = and o"'(0) = are 
satisfied. Let 

Sp = {a2,a3, ...,ap). (8) 

In order to obtain the local optimal values of its coefficients, we carry out the NM simplex method 
by using a simple initial value Sp = (0, ...,0, 50) for all p = 2,3, ...12, respectively. The result is 
summarized in Table 1. 

In general, from Table 1, it can be deduced that when p gets bigger, the number of iteration 
gets bigger. And \R\ improves when p increases from 2 to 10, but stopped improving when p > 10. 
The best result comes at p = 10 and 11 by giving |i?|=0.0031, only 24% of the optimal value 0.013 
obtained in [2]. The number of iteration are 1279 and 715, respectively. 

We also observe that: (a) in Table 1, some of the coefficients 03,04, ...,ap_i are not significant 
when compared to 02 and Op; (b) when a gets smaller more rapidly as r — s- and bigger more 
rapidly as r — > 1, the result gets much better. 



Table 1: Local optimal values of the coefficients of a defined by (6). Results are derived by the 
Nelder-Mead simplex method. At p = 6, 9, 12, the algorithm stops because the number of function 
evaluation has exceeded the preset limit of 2000. The actual number of iteration will be bigger if 
we increase the limit. 
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C. A Simpler Function Class for the PML Profile 

Based on the above observations, we try to define a simpler - in fact, better - function class for 
a: due to (a), let a/c = for /c = 3, — 1 to reduce complexity; due to (b), let the denominator 
of a in (6) be (1 — r), such that cr ^ as r —> and cr ^ oo as r ^ 1. After these changes, the 
shape of a satisfies the characteristic in (b) better. Then (6) becomes 

a2r2 + apTP 

a = —,p>2 (9 

1 — T 

where a{0) = and ct'(O) = still hold. 

We carry out the NM simplex method for a defined in (9) again with the initial value Sp = 
(0, 0, 50) for all p = 2,3, ...12, respectively. The local optimal values of its coefficients is listed in 
Table 2. 

From Table 2, we can see that after simplification, the maximal number of iteration is reduced 
to 157. (In Contrast to Table 1, the maximal number of iteration is more than 1370.) gets 
smaller when p increases from 2 to 8, but bigger when p > 9. The best result |i?|=0.0037 occurs 
at p = 8 and 9. It is bigger than |ii|=0.0031 obtained by (6), however, compared with its fast 
convergence, (9) is superior in efficiency and more practical than (6). In practice, we could choose 
a profile which has a lower order: for example, p = 5, the corresponding |i?|=0.0047. 

D. Comparison 

For more details, we plot as functions of 6 for four different a in Fig. 1 and Fig. 2. In Fig. 

1, we can see that a = (02''"^ + a8''"*)/(l — t) gives the lowest average of \R\ for 9 > 0.37r/2. In Fig. 

2, for very small angles {9 < 0.0057r/2), it again distinguishes itself from the other three by giving 
the lowest reflectivity. 

To illustrate that the NM method gives local minimizers, we plot as multivariable functions 
of a2 and ag using a = (a2r^ + a8T*)/(l — r) in Fig. 3. It can be observed that at 02 = 23.3 and 
as = 121.3, \R\ almost reaches its lowest value. 
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Fig. 1. Discrete reflectivity of four different PML profiles using the optimal values listed in Table 

1 and Table 2. For the dotted hne, S = 100.4. 
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Fig. 2. An enlargement of Fig. 1 for small angles. 



4. Conclusion 

For frequency domain simulations, a procedure for selecting the optimal PML profile is estab- 
lished based on minimizing the average discrete refiectivity. By extending the profile to a rational 
function class and adopting the Nelder-Mead simplex method to calculate the profile's coefficients, 
we reach a better numerical result. We also provide a simpler profile, which largely conserves the 
good properties of the former rational function. 

For further improvements, we may use a better function class for the PML profile, or improve 
the convergence property of the optimization method. In addition, we may study the impact of 
such proposed profile in a more practical example, for instance, optical wave propagating along 
an optical waveguide and hitting the PMLs in different angles, or how the guided and evanescent 
waves behave at a waveguide discontinuity, etc. 



Table 2: Local optimal values of the coefficients of a defined by (9). Results are derived by the 
Nelder-Mead simplex method. 
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Fig. 3. The average of reflectivity as functions of a2 and ag- The white arrow points at 02 = 23.3 

and as = 121.3. 
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